Core-softened potentials and the anomalous properties of water 
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We study the phase diagram of a system of spherical particles interacting in three dimensions 
through a potential consisting of a strict hard core plus a linear repulsive shoulder at larger distances. 
The phase diagram (obtained numerically, and analytically in a limiting case) shows anomalous 
properties that are similar to those observed in water. Specifically, we find maxima of density 
and isothermal compressibility as a function of temperature, melting with volume contraction, and 
multiple stable crystalline structures. If in addition a long range attraction between the particles 
is included, the usual liquid-gas coexistence curve with its critical point is obtained. But more 
interestingly, a first order line in the metastable fluid branch of the phase diagram appears, ending 
in a new critical point, as it was suggested to occur in water. In this way the model provides a 
comprehensive, consistent and unified picture of most of the anomalous thermodynamical properties 
of water, showing that all of them can be qualitatively explained by the existence of two competing 
equilibrium values for the interparticle distance. 
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I. INTRODUCTION 

Water is an anomalous substance in many respects.0 
Liquid water has a maximum as a function of temper- 
ature in both density and isothermal compressibility. It 
solidifies with volume increasing at low pressures, and the 
solid phase (ice) shows a remarkable variety of crystalline 
structures in different sectors of the pressure-temperature 
plane. Some of these properties are known from long 
ago, but their origin is still controversial. In an effort to 
rationalize these anomalous properties, the supercooled 
(metastable) sector of the phase diagram of liquid ^a- 
ter has received much attention in the last years nu It 
was observed that when appropriately cooled (using tech- 
niques for preventing crystallization) water becomes a 
viscous fluid with many properties (as heat capacity and 
isothermal compressibility) displaying a tendency that 
has suggested even .a thermodynamic singularity at some 
lower temperature.|j Although there is a limit of about 
235 K below which water cannot be cooled without crys- 
tallization, amorphous states of water at much lower 
temperatures can be obtained by different techniques. 
All these amorphous states are observed to correspond 
to one of two different structures (referred to as low- 
density amorphous -LDA- and high-density amorphous 
-HDA) that differ by about 20 per cent in density, which 
transform reversibly one into the other upon changes of 
pressure.B There is evidence that these amorphous states 
are thermodynamically connected with fluid water, al- 
though a direct verification is not possible due to recrys- 
tallization at intermediate temperatures!] 

The observation of LDA and HDA was an experimen- 
tal clue that led r to . th e proposal of the second critical 
point hypothesis .cft30 This hypothesis states that in 
the deeply supercooled region water can exist in two dif- 
ferent amorphous configurations, separated by a line of 
first order transitions. This line should end in a criti- 



cal point very much as the usual liquid-vapor line ends 
in a critical point. This hypothesis, in addition to obvi- 
ously explain the reversible transformation between LDA 
and HDA, provides a natural though phenomenological 
explanation for the anomalous behavior of density and 
isothermal compressibility. However, the very existence 
of the second critical point is known to be .not neccesary 
for the appearance of other anomalies ,E3Eil and the issue 
of what are the microscopic properties of water molecules 
that may produce the appearance of the second critical 
point are only poorly understood. In all cases it seems 
to be crucial the fact that water (because of the partic- 
ular form of its molecules and peculiarities of the hydro- 
gen bond) exhibits competition between more expanded 
structures (preferred at low pressures) and more com- 
pact ones (which are favored at high pressures). But it is 
not obvious to what extent this simple fact can be made 
responsible for all the anomalies of water, or if more sub- 
tle properties of the interactioiMiotential (in particular, 
cooperative hydrogen bonding)llM-a are crucial. 

Numerical simulations based on some of the avail- 
able pair potentials for the interaction between wa- 
ter molecules reproduce reasonably well many of its 
properties, although the systems that have been stud- 
ied are stroiis!^ limited in size, due to computational 
constraints OEJ These simulations only suggest the exis- 
tence of the second critical point, but up to now they were 
not able to prove its existence unambiguously. Other 
simplified and in some cases ad-hoc models have been 
devised to show the appearance of anomalous properties 
in the phase diagram. E3 Some of these models have a 
second critical point, but in these cases a global charac- 
terization of the phase diagram that includes all other 
anomalies has not bee achieved. In all cases, the models 
used have as a fundamental ingredient the competition 
between expanded, less dense structures and compressed, 
more dense ones. 
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It is the goal of the present work to show that a very 
simple model of spherical particles interacting trough a 
repulsive potential that possesses two different preferred 
equilibrium positions has a phase diagram in which a) 
lines of maxima for the density and the isothermal com- 
pressibility of the liquid exist; b) the fluid phase freezes 
with an increase in volume in some pressure range; c) 
the solid phase has multiple different crystalline struc- 
tures depending on P and T. When a long range van 
der Waals attraction is included on top of the previous, 
exclusively repulsive potential, the system d) preserves 
the anomalies existent in the non-attractive case; e) de- 
velops a liquid gas first order coexisting line that ends 
in a critical point in the usual fashion; f) depending on 
the strength of the attractive potential, a line of first or- 
der transitions separating two amorphous phases in the 
supercooled region appears. This line ends in a second 
critical point from which the line of maxima in isothermal 
compressibility starts. 

These statements will be justified mostly using numer- 
ical (Monte Carlo) techniques. However, for the deep su- 
percooled states, the long equilibration times make the 
numerical studies not completely reliable. In this case, 
the numerical results are supported by analytical calcu- 
lations in a limiting case of the interacting potential that 
shows neatly how the second critical point appears. 

The paper is organized as follows. In Section II we 
give a brief description of the interaction potential and 
the numerical technique. Section III focus on the dif- 
ferent stable crystalline configurations of the system. In 
Section IV we study the phase diagram of the fluid phase, 
both where it is thermodynamically stable and also in the 
supercooled region. Here we rely both in analytical calcu- 
lations and in simulations, and show how a second critical 
point can appear. In Section V we describe the melting 
of the most expanded solid structure and its anomalies. 
Finally, in Section VI we take all the results as a whole 
and comment upon their importance for the understand- 
ing of the properties of water. 

II. MODEL AND NUMERICAL DETAILS 

The interaction potential U(r) between particles that 
we will consider is chosen to be the hard-core plus lineac=. 
ramp potential originally studied by Stell and Hemmer.t£l 
The radius of the hard core is taken to be ro, and the 
ramp extends linearly from the value £o at r = ro up to 
at r = 7*1. In addition, a long range van der Waals 
attraction will be included through a global term in the 
energy per particle of the system proportional to — j/v , 
with v the specific volume and 7 a coefficient that rep- 
resents the total integrated strength of the attraction. 
The van der Waals term can be accounted for without 
its explicit inclusion in the simulations in the following 
way.Ej Since the free energy per particle contains the term 
Pv — j/v, when minimizing with respect to v the combi- 



nation P + 7/u 2 appears. So we will call P* = P + j/v 2 , 
and make the simulations in terms of P* , with 7 = 0. At 
the end, the self consistent replacement P* — > P + 7/u 2 
is made, and this provides the results for finite 7. 

Numerical simulations are performed at constant P, 
T, and N (the number of particles) by standard Monte 
Carlo techniques. Periodic boundary conditions are used 
in the three directions. The equilibrium volume at each 
pressure is reached by allowing the system size to increase 
or decrease through Monte Carlo movements that expand 
of contract all coordinates of the particles as well as the 
total size of the system. The rescaling is accepted of 
rejected depending on the energy change it involves. The 
contraction-expansion procedure is made independently 
for the three spatial coordinates, and the maximum ratio 
between the size of the system in different directions is 
limited to 1.2. 

A subtle but important technical change was intro- 
duced in the simulation procedure to be able to equi- 
librate the volume in the low temperature region. Let 
us think for instance of the T — case. If we increase 
P, as soon as two particles are at a distance ro from 
each other, the total volume gets stuck in the simula- 
tion because (due to the scheme adopted for doing vol- 
ume changes) the volume can be reduced further only 
if a global contraction of all coordinates reduces the en- 
ergy. But this contraction would bring the two particles 
(which arc already at a distance ro) at a distance lower 
than ro, then formally to a state of infinite energy, and 
so the trial movement is rejected. To avoid this problem 
we relax a bit the rigid hard core at ro , replacing it by a 
new linear ramp between r = and r — ro of the form 
U(r) = e [50(l — r/r ) + 1]. The last term is included 
in order to match smoothly the potential for r > ro. 
This modification was seen not to modify the behavior of 
the system, it just provides a convenient way of reaching 
the equilibrium values of the volume within a reasonable 
computing time in the low temperature regime. 

The property of the potential that renders it interest- 
ing for our problem is that depending on the external 
force acting on them, two particles prefer to be at dis- 
tance ro or r\ from each other, and the effect of this 
simple fact on the phase diagram is dramatic.Ej Already 
in the original papers about this potentialJLS it was real- 
ized that in fact the competition between configurations 
with particles at distances ro and T\ may produce the ap- 
pearance of polimorfism in the system. More precisely, 
in ID the system may exhibit many liquid phases when 
7^0, with sharp transitions between them. In 3D these 
transitions occur within the solid region of the phase di- 
agram, and it was suggested that they may still, hp .ob- 
servable as isostructural transitions of the solidJiflO In 
this connection we want to emphasize the following two 
points: i) isostructural transitions within the solid phase 
for this kind of potentials are usually preempted by the 
appearance of new intervening solid phases of different 
symmetry; ii) polimorfism of the liquid state as observed 
in the ID system also appears in 3D samples, but now 
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in the supercooled liquid state, and is the responsible for 
the existence of the second critical point. We address 
these two points in the next two sections. 

III. CRYSTALLINE CONFIGURATIONS AT T = 

Multiple crystalline structures for our potential arise 
from the competition between expanded and contracted 
structures. At T = the preferred configuration of 
the system will be the one that minimizes the enthalpy 
h = e + Pv. At low P the best way of minimizing h 
is first to minimize e, and then v. This is achieved in 
a close packed structure with first neighbors distance 
equal to r\. At very large P, in order to minimize h 
it is energetically more convenient to minimize v, and 
the structure becomes again a close packed one with first 
neighbors distance equal to 7*0. However, in the range in 
which these two configurations have approximately the 
same enthalpy, there are others which are more stable, 
having pairs of neighbor particles both at distances ro 
and 7"i . For our potential U (r) (and also for more gen- 
eral potentials) the existence of other stable-configura- 
tions can be demonstrated quite generally.E3c3 However, 
to tell safely which is the structure of lowest enthalpy 
for each P is a problem for which a closed solution is 
not known. The usual approach is to compare the en- 
thalpies of structures proposed beforehand, using simu- 
lated annealing techniques to guess the possible struc- 
tures. We refer to [22,23] for a detailed discussion of 
the two-dimensional case, and also for discussions on the 
neccesary conditions on the potential for different struc- 
tures to appear. 



FIG. 1. Crystalline structures of highest stability for dif- 
ferent values of P and 7, for n/ro = 1.75 and T = (adi- 
mensional values of P and 7 are defined as Po = ProejT 1 , 
and 70 = 7£^" 1 r ( 7 3 . The search was performed among struc- 
tures of the cubic (simple (labeled SC), face centered (FCC), 
and body centered), tetragonal (simple, and body centered), 
rhombohedral (RH), and hexagonal (simple (H), and close 
packed (HCP)) crystalline systems, with only one particle per 
unit cell (except for the HCP structure) . However, other more 
complex structures cannot be ruled out. 

Here we show only in Figure ^ the result of comparing 
(for 7*1/7*0 = 1.75) the enthalpies of particles arranged in 
Bravais lattices corresponding to cubic, tetragonal, rhom- 
bohedral, and hexagonal systems as a function of P and 
7. The structures searched are those that can be defined 
by no more than two parameters (that fix the form and 
size of the Bravais lattice). All of them were supposed 
to have only one particle per unit cell, except for HCP 
(which has two), that was included due to its known sta- 
bility. We find five different crystalline configurations as 
a function of P. Note also how the increasing of the 
van der Waals attraction moves all the borders between 
structures to lower pressures, since those that are stable 
at higher P are always more compact, and thus become 
more stable in the presence of the van der Waals term. 
It must be kept in mind that there may be other config- 
urations (corresponding to other crystalline systems, or 
with more complex unit cells) with lower enthalpy. In 
fact, they are likely to occur, as for instance in the 2D 
case crystalline structures with up to five particles per 
unit cell appear.E3o Only a thorough numerical work 
can determine all possible structures. 
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IV. PROPERTIES IN THE FLUID AND 
SUPERCOOLED REGIONS 

In the previous section we saw that at T = there is 
a sequence of solid phases interpolating between the low- 
est density and highest density ones. At each transition 
(as P is increased) a finite fraction of particles that were 
at distance r\ from each other passes to be at distance 
ro. Each of these rearrangements involves a change of 
symmetry, and thus the appearance of a new crystalline 
structure. The picture is different in the metastable, dis- 
ordered sheet of the phase diagram. At very low pressures 
the particles behave as hard spheres with radius r\j1. 
Hard spheres are known to have a maximum density of 
random packing, corresponding to a volume per particle 
v\ than in our case is V\ ~ 0.808^. Being the dens- 
est disordered structure possible for hard spheres, this is 
also the thermodynamically stable amorphous configura- 
tion of the system when T = 0, namely, the one which 
minimizes the enthalpy. When P — > 00 the linear ramp 
of U(r) is irrelevant to calculate the free energy, and the 
thermodynamically stable configuration is again a ran- 
dom packing of spheres, now with radius ro/2. As in the 
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crystalline case, the nearest neighbor distance between 
particles must collapse from n to rg as a function of P. 
The crucial question is whether this collapse is discontin- 
uous at some well defined P (or even if there are more 
than one transitions at different values of P) or if it is just 
a smooth crossover. We address the issue in the follow- 
ing two subsections, first analytically (when ri/ro —> oo) 
and then numerically, but let us quote briefly the answer 
to this question in advance. For the case of no van der 
Walls attraction the behavior of the specific volume v as 
a function of P is smooth at any finite temperature. But 
there is a range of pressures in which -^p is anomalously 
large. This fact is enough for the van der Waals interac- 
tion to produce (if sufficiently strong) the appearance of 
a metastable critical point and a line of first order tran- 
sitions between two disordered structures with a finite 
difference in density. 



A. The limit ri/ro — > oo 

We will start by considering only the repulsive part of 
the potential (i.e., 7 = 0). As we already said, at very 
low pressures the particles behave as hard spheres with 
radius ri/2. At T = the enthalpy per particle of this 
configuration is h — Pv\ . This is the thermodynamically 
stable state upon increasing P up to the point where it 
is energetically more convenient to overlap neighbor par- 
ticles in pairs. The structure will now be similar to the 
low pressure one, but with two particles overlapped on 
each position (see a sketch of this fact in Figure 0) . The 
enthalpy of this configuration is h' = Pvi/2 + £o/2, since 
now the total volume of the system is reduced in a fac- 
tor 2, and an energy eo must be counted for each pair 
of particles. The pressure at which h = h! determines 
the transition pressure Ptr. — £o/ v i- Close to the point 
(Ptr, T = 0) of the phase diagram, we can calculate ap- 
proximately the free energy of the system in the following 
way. Let us suppose we have N particles, n of them in 
non overlapped positions and n' pairs of overlapped par- 
ticles (N — n + 2n'). The configurational free energy of 
the system may be written as (higher than double over- 
laps that will ocurr at higher pressures are dismissed) 



F = [Pv' - Ts Hb (v')](n + n') + e n' - Ts 1 



(1) 



where v' = V/ (n + n 1 ) [note that v' is not the specific vol- 
ume, which instead is given by v = V/(n + 2n')], s HS (v 1 ) 
is the entropy per particle of hard spheres with radius 
n on the metastable sheet, and s 1 is the configurational 
entropy for choosing which particles will be in pairs, and 

n + n' 



which ones will be singled 
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ing v' and n' as independent variables for minimizing F 
we obtain the equations 



P _ d. 
T ~ ~ 



,HS t~i 



(«') 



dv' 



(2) 



p v' H-9/-A £o , , (l-2n'/Nr 

(we use v' in this case with 7 = 0, to distinguish from 
the 7^0 case). 

The first one is the equation of state of hard spheres 
in the metastable region. We will use for it the following 
expression provided by SpeedyEa 



P _ 2.65fc B 
T v' — vi 



(4) 



For given values of P and T, v' is determined from this 
equationj-and the value obtained is used to determine nl 
from (g)Ej 
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FIG. 2. Sketch of two dimensional configurations of parti- 
cles in a glassy state at T = for ri/ro large, although finite, 
below and above the transition pressure Ppr. Drawings were 
made so as to emphasize that for P > Ptr (b) the structure 
is similar to that at P < Ptr (a) but with two particles per 
site, instead of one (this fact is strictly valid if ri/ro — > 00). 
The density of the system is roughly twice in (b) than in (a). 
The energy per particle is in (a) and Eo/2 in (b). 

The volume per particle of the system is v = V/N — 
v'(n + n')/N. Although v' has a behavior on P and T 
that is the same as for hard spheres, the (n + n') factor 
(that takes the value N when T = 0, P < P T r, and N/2 
when T — 0, P > Ptr.) makes the behavior of v on T 
and P be non trivial. We see the surface v(P,T) that is 
obtained from equations (§), (g), and (@) in Fig. || At 
T = we obtain the expected result, with v(P,T = 0) 
passing from v\ to v\/2 at the transition pressure Ptr- 
But at any finite T the entropy transforms this jump in 
a smooth crossover. The point P = Ptr, T = is for 
this system the metastable critical point. 

The inclusion of a finite long range attraction through 
a non zero 7 produces the critical point to move into 
the T > region. The mechanism is identical to the 
one that produces the appearance of the usual liquid-gas 
coexistence curve. We must replace P by P + 7/f 2 in 
expressions (H), (g), and (g) to take account of the van 
der Waals attraction. A singularity in v(P) at a finite 
temperature exists if dv/dP becomes negative (this is 
the signature of a van der Waals loop, and thus of a first 
order transition). Since 



v(P) =i(P + 1 /v 2 ), 



(5) 
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we can calculate dv / dP as 



dv/dP 



dv{x)/dx 



1 



27 dv(x) 



(6) 



X=P-\-~f/ V 2 



We see that a singularity occurs if dv(P)/dP is larger 
(in absolute value) than |-. This always happens in our 
model close enough to Ptr, T = 0. In Figure |I] the func- 
tion v(P, T) , calculated using the self consistency condi- 
tion (||), with 71 = 0.1 (71 = 7£ ~ 1 r 1 ~ 3 )is shown. 



We see the location in the P-T diagram of the extrema 
of v and Kt as a function of temperature in Figure |^. 
We also see in this figure the first order line that appears 
due to the van der Waals attraction, ending in the crit- 
ical point C, and also the two spinodal lines that mark 
the limit of metastability of the two phases on both sides 
of the first order line. Note that the singularity of the 
thermodynamic properties at the critical point manifests 
itself in anomalous properties (of v and Kt) that can be 
detected at higher temperatures. 




FIG. 3. The surface v(P,T) for 71 = (71 = ye^r' 3 ), 
when n/ro — > 00. Only the T — isotherm is singular, at 
P = Ptr (T is measured in units of kg so, Pi = Pr^e^ 1 ). 




T 



FIG. 4. Same as Figure (H), but with 71 = 0.1. Now a 
discontinuity exists for all T^0.37. 

The rapid change in v as a function of P close to the 
critical point is the responsible for the anomalous behav- 
ior of v and the isothermal compressibility Kt = — ~ §p ■ 




FIG. 5. The locus of extrema of v and Kt, calculated from 
the v(P,T) function, for 71 = (upper curves) and 71 = 0.1 
(lower curves). For finite 71 a first order line ending in a 
critical point C' appears. The dashed-dotted lines are the 
spinodals of this first order transition. 

The analytical treatment of the case Tx/tq — > 00 pro- 
vides insight into the appearance of the second critical 
point in the phase diagram of water. In fact, the anoma- 
lies in v and Kt exist even for exclusively repulsive po- 
tentials. It is the van der Walls attraction that brings 
the critical point to a finite temperature, in the same 
way that it is this attraction (or a more realistic finite 
range one) that generates the familiar liquid-gas coexis- 
tence line. Now we will see how much of this scenario 
remains for finite ri/fo- 

B. Numerical results for finite ri/ro 

When ri/ro is finite, no analytical calculation seems to 
be possible to tell the existence or not of the metastable 
critical point. But guided by the previous findings, we 
can more safely interpret the numerical results. 

In Figure |^ we show the results of numerical simula- 
tions for n/ro = 1.5. Rapid runs (2000 steps per temper- 
ature) were made by decreasing the temperature at dif- 
ferent values of P in a system of 197 particles. The rapid 
cooling allows to reach the supercooled states without 



5 



crystallization except in the continuous-dashed region.Ea 
The curves shown are averages over 20 different runs. 
Only the points in which the system did not crystallize 
and displays well reproducible values for the density are 
shown. In spite of this, since the runs were rapid and the 
low temperature states are highly viscous, we can rise 
some doubts about the final state reached for T — > 0. It 
might be that we are observing some frozen configura- 
tion typical of larger T. To answer this point we made 
runs at a low temperature (T = 0.01) increasing and de- 
creasing pressure (Fig. [?]) . The results of this simulation 
show evident effects of hysteresis due to the glassiness of 
the states. This hysteresis was seen not to be greatly re- 
duced by decreasing the rate of temperature change in a 
factor ten. But anyway the hysteresis path encloses the 
values of v obtained by decreasing T at fixed P (large 
symbols in Fig. (jj), and we have also checked that the ra- 
dial distribution functions are comparable in both cases. 
The finding of essentially the same results when we ar- 
rive from different paths in the P-T plane is an indication 
that in fact these are thermodynamic values. 
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T 0.20 



FIG. 6. Specific volume as a function of T for different 
values of Pq (Pq — PtqEq 1 ) from 0.2 to 2.8 in steps of 0.2, from 
the numerical simulations with n/ro = 1.5. Simulations were 
done reducing rapidly the temperature from the fluid state. 
Averages over 20 runs are shown. The continuous-dashed 
region corresponds to points where the system crystallization 
could not be avoided even in these rapid runs, and are thus 
not included. Within the region limited by the dashed line 
dv/dT is negative. 



are sufficient for the van der Walls attraction to induce 
the appearance of a critical point, if 7 is large enough. 
In fact, we find that for 70 = 7p r — 4.2 a critical point 
enters the phase diagram at T = 0, Pa ~ 0.65. The 
location of the critical point as a function of 7 can be 
determined from data as those of Fig. ||by requiring that 
dv/dP — > 00, and d 2 v/dP 2 — > 00 calculated according 
to Eq. (||). We find in our case that for 70 slightly larger 
than 7Q r the critical point position can be estimated as 



T cr ~ 0.07(7o - 4.2) 
P cr ~ 0.65 -0.35(7o -4.2). 



(7) 
(8) 




4.0 Po 5.0 



FIG. 7. Specific volume as a function of P at T = from 
the simulations shown in Fig. ^ (large symbols) and from 
a single run at T = 0.01 (20000 steps per temperature) in- 
creasing and decreasing P (small symbols. The path of the 
simulation is indicated by the arrows). The dashed region 
contains no large symbols since here crystallization could not 
be avoided when reducing T. Note that the limiting values of 
v/vq at very high and low P are respectively 0.808 and 2.727. 

The locus of the anomalies of v and Kt also move with 
7. We note that if 7 is such that the critical point exists, 
the line of Kt maxima necessarily ends at the critical 
point (since at the critical point Kt — > 00). For the ex- 
trema of v this is not necessarily so, although it is known 
that the_anomalies in Kt and v are thermodynamically 
related.B 



There is no sign in Fig. (in contrast to the fi/ro — > 00 
case) of an abrupt jump in v as a function of P at T = 0, 
all that remains is a value of pressure with a maximum in 
dv{P,T = 0)/dP (as is seen in Fig. g, close to P = 2). 
Around this value Kt has a maximum (as a function 
of T) at T = 0, whereas v has maxima and minima at 
finite temperature, as can be seen from Fig. ^| in the 
range 1.0 < Pq < 1.8 (on the dashed line). These facts 



V. CHARACTERISTICS OF MELTING 

In this section we show results of numerical simula- 
tions that focus on the melting _af the most expanded 
of the solid phases of our systemcj (which is the equiv- 
alent of ice Ih in water). Figure @ shows the specific 
volume v as a function of T for different values of P ob- 
tained in slow simulations decreasing and increasing T in 
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a system of 216 particles. The hysteresis upon heating 
and cooling embraces the position of the thermodynamic 
melting transition temperature. In all the range of P 
indicated in this figure, the system freezes into one and 
the same solid configuration, corresponding to a dense 
stacking of triangular planes, in which each particle has 
twelve nearest neighbors at distance n (the dispersion 
in the limiting value of v when T — > is due to a few 
defects that remain in the solid structure). Upon increas- 
ing T, v increases for Pq < 0.95 (which is of course the 
standard behavior), but decreases for larger Pq. This de- 
crease is driven by the possibility of particles of being at 
distances smaller than r% from each other. Depending 
on P, the tendency of particles to become closer (gain- 
ing energy from the Pv term) may be higher than the 
entropic tendency to increase v. In the same way, at the 
lowest pressures, the solid melts by increasing its volume, 
whereas at the largest pressures shown in Fig. || it melts 
by reducing its volume. This is consistent with the form 
of the solid-liquid border in the P-T plane that is seen 
in Figure ^|, which has positive derivative at low P, but 
has negative derivative at larger P. 




0.3 j 



FIG. 9. The melting line of the lowest pressure solid struc- 
ture, for 70 = 0, 3, 5, and 6. There is also a liquid-gas coexis- 
tence line for all 70 7^ that in the scale of the figure cannot 
be distinguished from the Po — axis. This line ends in the 
critical point indicated for each case by a square. 
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FIG. 8. The specific volume v as a function of T for dif- 
ferent values of Po, from 0.1 to 1.15, in steps of 0.05 for a 
system of 216 particles. Each curve was vertically displaced 
by a term — 2Pq to allow a better visualization. Hysteresis is 
the result of successive cooling and heating. Note the normal 
melting at low Po, and the anomalous one for Po>0.95. 



In the same Figure y we see the modification of the 
phase diagram when we consider the van der Waals at- 
traction. Any finite value of 7 makes a liquid-gas first 
order line appearJlS In the scale of Fig. |j this critical 
line cannot be distinguished from the P — axis, only 
the critical temperature is indicated. In addition, the 
whole solid-fluid coexistence line basically moves down 
with 7q. If 70 5.5 the triple point that is defined is 
'standard', in the sense that the slope of the solid-liquid 
coexistence line is positive at the triple point. For larger 
70 the triple point is 'anomalous' (the slope of the solid- 
liquid coexistence line is negative). 



VI. SUMMARY AND CONCLUSIONS 

We studied the phase diagram of a model of spherical 
particles with pairwise interactions, consistent of a hard 
core at a distance r$ plus a repulsive linear shoulder that 
extends up to distance r\. This potential favors the par- 
ticles to be in one or the other (depending on P) of the 
two different equilibrium distances and r\. On top 
of that, a long range van der Waals attraction was also 
included. 

The solid phase of the system exhibits polimorfism. 
Namely, there are different sectors of the P-T phase di- 
agram in which the crystalline structure of the system is 
different. This behavior is observed even in the case of 
non attractive part in the potential (i.e., 7 = 0). 

The fluid phase part of the phase diagram has the fol- 
lowing characteristics. At low pressures particles prefer 
to be at distances ~ r x from each other, whereas at high 
pressures the typical distance is r$ (< r%). This implies a 
crossover region of P with an anomalously large isother- 
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mal compressibility Kt- When we include the van der 
Walls attraction (7 7^ 0) the anomaly in Kt may become 
(if 7 is sufficiently large) a first order transition line (sim- 
ilar to the liquid-gas coexistence line). This first order 
line starts from a finite P at T = 0, and ends in a critical 
point at finite T and P. From this point the line of Kt 
maxima continues towards larger values of T. There are 
also anomalies in the density of the system, which has 
extrema in a locus that, whereas it does not necessar- 
ily touch the critical point, appears in the region that is 
influenced by the existence of it. 

For 7 = 0, the melting line of the most expanded solid 
structure in the P-T plane has positive derivative at low 
pressures, but is reentrant at higher P. This reentrant 
behavior is associated (through the Clausius-Clapeyron 
equation) to a melting with density increasing. When 
the van der Waals attraction is included, a liquid-gas first 
order line appears, that ends in a critical point as usual. 
This liquid-gas line defines a triple point where it meets 
the fluid-solid line. For small 7 the slope of the solid- 
liquid line at the triple point is positive, but it becomes 
negative if 7 is large enough. 

Our model, although very simple, has many of the 
properties that characterize water as an anomalous fluid, 
and gives insight into the properties of real water. Ac- 
tually, the simplicity of the model allows to single out 
the crucial characteristic that produces all the anomalies, 
without the complications introduced by non-spherical 
interactions and cooperative hydrogen bonding in real 
water. This characteristic is the existence in the inter- 
atomic potential of two different equilibrium distances 
for the particles. From all our results it is difficult to 
elude the claim that there must be an effective descrip- 
tion of the interaction in water in which two diggerent 
distances compete as being the most stable one. In fact, 
it would be even more daring to say that all the similar- 
ities we found are accide ntal. Although there is evidence 
favoring this view,ljO'E3 the complications added by the 
peculiarities of water molecules has made this point very 
disputed.E3 

Among all anomalous properties of water, the existence 
of the second critical point is the one that is not fully 
proven to ocurr, and also the one that has been most 
elusive to adress numerically in previous studies. Our 
model shows that its existence is a consequence of the 
effect of the attractive part of the potential on a system 
that (due to peculiarities of the interaction) possesses 
anomalously large values of Kt at some pressure. At 
this point experimental evidence about the existence of 
two amorphous phases (LDA and HDA) that transform 
reversibly into each other seems crucial to indicate that in 
water, the attraction between molecules is strong enough 
to bring the second critical point into existence. 

From a more fundamental point of view, we note that 
our model has essentially two free parameters, the ratio 
between equilibrium distances r\/ra and the strength of 
the van der Waals attraction 7. Other characteristics, 
as if the ramp between rg and r\ is linear or not, are 



only marginal for the phase diagram that is obtained.EHI 
An important result of our study is the fact that these 
two parameters determine the phase behavior of the fluid 
phase both in the zone where the fluid is stable, and also 
in the deeply supercooled region. If water admits a siixu 
ilar effective representation in term of two parameters,E3 
then these could be extracted from fitting experimental 
data in the high temperature region, and then used, rely- 
ing on our model, to predict the supercooled part of the 
phase diagram, in particular the existence and location 
of the second critical point. This means that the present 
model may even be of quantitative importance. Work on 
this direction is under way. 
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